## Callin Switzer
## 17 Nov 2017
## Multilevel model to visualize bees'
## behavior on the artificial pollen system
#install packages
ipak <- function(pkg){
new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}
packages <- c("ggplot2", "reshape2", 'lme4', 'sjPlot', "multcomp", "plyr", "effects", "viridis")
ipak(packages)
## Loading required package: ggplot2
## Loading required package: reshape2
## Loading required package: lme4
## Loading required package: Matrix
## Loading required package: sjPlot
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.2.10
## Current Matrix version is 1.2.11
## Please re-install 'TMB' from source or restore original 'Matrix' package
## Loading required package: multcomp
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
##
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
##
## geyser
## Loading required package: plyr
## Loading required package: effects
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
## Loading required package: viridis
## Loading required package: viridisLite
## ggplot2 reshape2 lme4 sjPlot multcomp plyr effects viridis
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
# set ggplot theme
theme_set(theme_classic() + theme(axis.text=element_text(colour="black")))
# define data and figure directories
dataDir <- "/Users/cswitzer/Dropbox/SonicationBehavior/SonBehData"
figDir <- "/Users/cswitzer/Dropbox/SonicationBehavior/SonBehFigs"
print(paste("last run ", Sys.time()))
## Warning in as.POSIXlt.POSIXct(x, tz): unknown timezone 'zone/tz/2017c.1.0/
## zoneinfo/America/Los_Angeles'
## [1] "last run 2018-01-11 17:39:46"
print(R.version)
## _
## platform x86_64-apple-darwin15.6.0
## arch x86_64
## os darwin15.6.0
## system x86_64, darwin15.6.0
## status
## major 3
## minor 4.2
## year 2017
## month 09
## day 28
## svn rev 73368
## language R
## version.string R version 3.4.2 (2017-09-28)
## nickname Short Summer
Read in data and double check it
sl <- read.csv(file.path(dataDir, '01_CombinedTrials_cleaned.csv'))
# make sure hive is a factor
sl$hive = as.factor(sl$hive)
unique(sl$hive)
## [1] 4 3 5
## Levels: 3 4 5
# number of trials per bee
xtabs(~sl$beeCol)
## sl$beeCol
## blue gold goldred green
## 36 54 3 13
## lime limeblue limegold limegreen
## 28 530 54 50
## limeorange limepink limepurple limepurpleyellow
## 84 51 352 58
## limered limesilver limewhite limeyellow
## 3474 3 101 52
## orange orangeblue orangegreen orangepink
## 501 50 50 50
## orangepurple purple redblue redgreen
## 50 9 673 530
## redpink redpurple silver white
## 677 218 26 283
## whiteblue whitegold whitegreen whiteorange
## 2786 1225 39 1122
## whitepink whitepurple whitered whiteyellow
## 1184 55 1344 157
## yellowblue yellowgreen yelloworange yellowpink
## 1587 532 2113 2309
## yellowpurple yellowred
## 1243 547
hist(xtabs(~sl$beeCol))
max(xtabs(~sl$beeCol))
## [1] 3474
mean(xtabs(~sl$beeCol))
## [1] 578.6429
sd(xtabs(~sl$beeCol))
## [1] 829.7747
# number of visits for first trial
xtabs(~sl$beeCol[sl$trialNum == 1])
## sl$beeCol[sl$trialNum == 1]
## blue gold goldred green
## 36 54 3 9
## lime limeblue limegold limegreen
## 28 50 54 50
## limeorange limepink limepurple limepurpleyellow
## 33 51 53 58
## limered limesilver limewhite limeyellow
## 50 3 50 52
## orange orangeblue orangegreen orangepink
## 85 50 50 50
## orangepurple purple redblue redgreen
## 50 9 49 51
## redpink redpurple silver white
## 56 55 26 48
## whiteblue whitegold whitegreen whiteorange
## 56 48 39 54
## whitepink whitepurple whitered whiteyellow
## 64 55 59 55
## yellowblue yellowgreen yelloworange yellowpink
## 54 51 54 58
## yellowpurple yellowred
## 54 57
hist(xtabs(~sl$beeCol[sl$trialNum == 1]))
mean(xtabs(~sl$beeCol[sl$trialNum == 1]))
## [1] 46.92857
sd(xtabs(~sl$beeCol[sl$trialNum == 1]))
## [1] 16.37725
# number of visits for second trial
xtabs(~sl$beeCol[sl$trialNum == 2])
## sl$beeCol[sl$trialNum == 2]
## blue gold goldred green
## 0 0 0 4
## lime limeblue limegold limegreen
## 0 54 0 0
## limeorange limepink limepurple limepurpleyellow
## 51 0 243 0
## limered limesilver limewhite limeyellow
## 281 0 51 0
## orange orangeblue orangegreen orangepink
## 0 0 0 0
## orangepurple purple redblue redgreen
## 0 0 60 53
## redpink redpurple silver white
## 58 90 0 67
## whiteblue whitegold whitegreen whiteorange
## 314 66 0 260
## whitepink whitepurple whitered whiteyellow
## 143 0 195 47
## yellowblue yellowgreen yelloworange yellowpink
## 128 50 176 142
## yellowpurple yellowred
## 346 56
# number of trials per bee and treatment
xtabs(~sl$beeCol + sl$trt)
## sl$trt
## sl$beeCol full high low
## blue 36 0 0
## gold 54 0 0
## goldred 3 0 0
## green 13 0 0
## lime 28 0 0
## limeblue 530 0 0
## limegold 54 0 0
## limegreen 50 0 0
## limeorange 84 0 0
## limepink 51 0 0
## limepurple 109 0 243
## limepurpleyellow 58 0 0
## limered 50 0 3424
## limesilver 3 0 0
## limewhite 101 0 0
## limeyellow 52 0 0
## orange 85 416 0
## orangeblue 50 0 0
## orangegreen 50 0 0
## orangepink 50 0 0
## orangepurple 50 0 0
## purple 9 0 0
## redblue 673 0 0
## redgreen 530 0 0
## redpink 56 0 621
## redpurple 55 163 0
## silver 26 0 0
## white 283 0 0
## whiteblue 56 0 2730
## whitegold 48 0 1177
## whitegreen 39 0 0
## whiteorange 54 0 1068
## whitepink 135 1049 0
## whitepurple 55 0 0
## whitered 59 1285 0
## whiteyellow 157 0 0
## yellowblue 54 1533 0
## yellowgreen 532 0 0
## yelloworange 54 2059 0
## yellowpink 58 0 2251
## yellowpurple 54 0 1189
## yellowred 547 0 0
colSums(xtabs(~sl$beeCol + sl$trt))
## full high low
## 5095 6505 12703
xtabs(~sl$beeCol + sl$IT_imputed)
## sl$IT_imputed
## sl$beeCol 3.01 3.4 3.48 3.61 3.67 3.74 3.75 3.76 3.82 3.87 3.96
## blue 0 0 0 36 0 0 0 0 0 0 0
## gold 0 0 0 0 0 0 0 0 0 0 0
## goldred 3 0 0 0 0 0 0 0 0 0 0
## green 0 0 0 0 0 0 0 0 0 0 0
## lime 0 0 0 0 0 0 0 0 0 0 0
## limeblue 0 0 0 0 0 0 0 0 0 0 0
## limegold 0 0 0 0 0 0 0 0 0 0 0
## limegreen 0 0 0 0 0 50 0 0 0 0 0
## limeorange 0 0 0 0 0 0 0 0 0 0 0
## limepink 0 0 0 0 0 0 0 0 0 0 0
## limepurple 0 0 352 0 0 0 0 0 0 0 0
## limepurpleyellow 0 0 0 0 0 0 58 0 0 0 0
## limered 0 0 0 0 0 0 0 0 0 0 0
## limesilver 0 3 0 0 0 0 0 0 0 0 0
## limewhite 0 0 0 0 0 0 0 101 0 0 0
## limeyellow 0 0 0 0 52 0 0 0 0 0 0
## orange 0 0 0 0 0 0 0 0 0 0 0
## orangeblue 0 0 0 0 0 0 0 0 0 0 0
## orangegreen 0 0 0 0 0 0 0 0 0 0 0
## orangepink 0 0 0 0 0 0 0 0 0 0 50
## orangepurple 0 0 0 0 0 0 0 0 0 0 0
## purple 0 0 0 0 0 0 0 0 0 0 0
## redblue 0 0 0 0 0 0 673 0 0 0 0
## redgreen 0 0 0 0 0 0 0 0 0 0 0
## redpink 0 0 0 0 0 0 0 0 0 0 0
## redpurple 0 0 0 0 0 0 0 0 0 0 0
## silver 0 0 0 0 0 0 0 0 0 0 0
## white 0 0 0 0 0 0 0 0 0 0 0
## whiteblue 0 0 0 0 0 0 0 0 0 0 0
## whitegold 0 0 0 0 0 0 0 0 0 0 0
## whitegreen 0 0 0 0 0 0 0 0 0 0 0
## whiteorange 0 0 0 0 0 0 0 0 0 0 0
## whitepink 0 0 0 0 0 0 0 0 0 0 0
## whitepurple 0 0 0 0 0 0 0 0 0 0 0
## whitered 0 0 0 0 0 0 0 0 0 0 0
## whiteyellow 0 0 0 0 0 0 0 0 0 0 0
## yellowblue 0 0 0 0 0 0 0 0 0 0 0
## yellowgreen 0 0 0 0 0 0 0 0 0 0 0
## yelloworange 0 0 0 0 0 0 0 0 0 0 0
## yellowpink 0 0 0 0 0 0 0 0 0 2309 0
## yellowpurple 0 0 0 0 0 0 0 0 1243 0 0
## yellowred 0 0 0 0 0 0 0 0 0 0 0
## sl$IT_imputed
## sl$beeCol 3.98 3.99 4 4.02 4.03 4.05 4.07 4.12 4.13 4.14 4.15
## blue 0 0 0 0 0 0 0 0 0 0 0
## gold 0 0 0 0 0 54 0 0 0 0 0
## goldred 0 0 0 0 0 0 0 0 0 0 0
## green 0 0 0 0 0 0 0 0 0 0 0
## lime 0 0 0 28 0 0 0 0 0 0 0
## limeblue 0 0 0 0 0 0 0 0 0 0 0
## limegold 0 0 0 0 0 0 0 0 0 0 0
## limegreen 0 0 0 0 0 0 0 0 0 0 0
## limeorange 84 0 0 0 0 0 0 0 0 0 0
## limepink 0 0 0 0 0 0 0 0 51 0 0
## limepurple 0 0 0 0 0 0 0 0 0 0 0
## limepurpleyellow 0 0 0 0 0 0 0 0 0 0 0
## limered 0 0 0 0 0 0 0 0 0 0 3474
## limesilver 0 0 0 0 0 0 0 0 0 0 0
## limewhite 0 0 0 0 0 0 0 0 0 0 0
## limeyellow 0 0 0 0 0 0 0 0 0 0 0
## orange 0 501 0 0 0 0 0 0 0 0 0
## orangeblue 0 0 0 0 0 0 0 0 0 0 0
## orangegreen 0 0 0 0 0 0 0 0 0 0 0
## orangepink 0 0 0 0 0 0 0 0 0 0 0
## orangepurple 0 0 0 50 0 0 0 0 0 0 0
## purple 0 9 0 0 0 0 0 0 0 0 0
## redblue 0 0 0 0 0 0 0 0 0 0 0
## redgreen 0 0 0 0 0 0 0 530 0 0 0
## redpink 0 0 0 0 0 0 0 0 0 0 0
## redpurple 0 0 0 0 0 0 0 0 0 0 0
## silver 0 0 0 0 0 0 0 0 0 0 0
## white 0 0 0 0 283 0 0 0 0 0 0
## whiteblue 0 0 0 0 0 0 0 0 0 0 0
## whitegold 0 0 0 0 0 0 0 0 0 0 0
## whitegreen 0 0 39 0 0 0 0 0 0 0 0
## whiteorange 0 0 0 0 0 0 0 0 0 0 0
## whitepink 0 0 0 0 0 0 0 0 0 0 0
## whitepurple 0 0 0 0 0 0 0 0 0 0 0
## whitered 0 0 0 0 0 0 0 0 0 0 0
## whiteyellow 0 0 0 0 0 0 0 0 0 0 0
## yellowblue 0 0 0 0 0 0 0 0 0 0 0
## yellowgreen 0 0 0 0 0 0 0 0 0 532 0
## yelloworange 0 0 0 0 0 0 0 0 0 0 0
## yellowpink 0 0 0 0 0 0 0 0 0 0 0
## yellowpurple 0 0 0 0 0 0 0 0 0 0 0
## yellowred 0 0 0 0 0 0 547 0 0 0 0
## sl$IT_imputed
## sl$beeCol 4.21 4.23 4.24 4.31 4.34 4.37 4.42 4.45 4.46 4.59 4.61
## blue 0 0 0 0 0 0 0 0 0 0 0
## gold 0 0 0 0 0 0 0 0 0 0 0
## goldred 0 0 0 0 0 0 0 0 0 0 0
## green 13 0 0 0 0 0 0 0 0 0 0
## lime 0 0 0 0 0 0 0 0 0 0 0
## limeblue 0 0 0 0 0 0 0 0 0 0 0
## limegold 0 0 0 54 0 0 0 0 0 0 0
## limegreen 0 0 0 0 0 0 0 0 0 0 0
## limeorange 0 0 0 0 0 0 0 0 0 0 0
## limepink 0 0 0 0 0 0 0 0 0 0 0
## limepurple 0 0 0 0 0 0 0 0 0 0 0
## limepurpleyellow 0 0 0 0 0 0 0 0 0 0 0
## limered 0 0 0 0 0 0 0 0 0 0 0
## limesilver 0 0 0 0 0 0 0 0 0 0 0
## limewhite 0 0 0 0 0 0 0 0 0 0 0
## limeyellow 0 0 0 0 0 0 0 0 0 0 0
## orange 0 0 0 0 0 0 0 0 0 0 0
## orangeblue 0 0 0 0 0 0 0 0 0 0 50
## orangegreen 0 0 0 0 0 50 0 0 0 0 0
## orangepink 0 0 0 0 0 0 0 0 0 0 0
## orangepurple 0 0 0 0 0 0 0 0 0 0 0
## purple 0 0 0 0 0 0 0 0 0 0 0
## redblue 0 0 0 0 0 0 0 0 0 0 0
## redgreen 0 0 0 0 0 0 0 0 0 0 0
## redpink 677 0 0 0 0 0 0 0 0 0 0
## redpurple 0 0 0 0 218 0 0 0 0 0 0
## silver 0 26 0 0 0 0 0 0 0 0 0
## white 0 0 0 0 0 0 0 0 0 0 0
## whiteblue 0 0 0 0 2786 0 0 0 0 0 0
## whitegold 0 0 0 0 0 0 0 0 0 0 0
## whitegreen 0 0 0 0 0 0 0 0 0 0 0
## whiteorange 0 0 0 0 0 0 0 0 1122 0 0
## whitepink 0 0 0 0 0 0 0 0 0 1184 0
## whitepurple 0 0 0 0 0 0 0 55 0 0 0
## whitered 0 0 0 0 0 0 1344 0 0 0 0
## whiteyellow 0 0 157 0 0 0 0 0 0 0 0
## yellowblue 0 0 0 0 0 0 0 0 0 0 0
## yellowgreen 0 0 0 0 0 0 0 0 0 0 0
## yelloworange 0 0 0 0 2113 0 0 0 0 0 0
## yellowpink 0 0 0 0 0 0 0 0 0 0 0
## yellowpurple 0 0 0 0 0 0 0 0 0 0 0
## yellowred 0 0 0 0 0 0 0 0 0 0 0
## sl$IT_imputed
## sl$beeCol 4.63 4.69 4.94
## blue 0 0 0
## gold 0 0 0
## goldred 0 0 0
## green 0 0 0
## lime 0 0 0
## limeblue 0 530 0
## limegold 0 0 0
## limegreen 0 0 0
## limeorange 0 0 0
## limepink 0 0 0
## limepurple 0 0 0
## limepurpleyellow 0 0 0
## limered 0 0 0
## limesilver 0 0 0
## limewhite 0 0 0
## limeyellow 0 0 0
## orange 0 0 0
## orangeblue 0 0 0
## orangegreen 0 0 0
## orangepink 0 0 0
## orangepurple 0 0 0
## purple 0 0 0
## redblue 0 0 0
## redgreen 0 0 0
## redpink 0 0 0
## redpurple 0 0 0
## silver 0 0 0
## white 0 0 0
## whiteblue 0 0 0
## whitegold 1225 0 0
## whitegreen 0 0 0
## whiteorange 0 0 0
## whitepink 0 0 0
## whitepurple 0 0 0
## whitered 0 0 0
## whiteyellow 0 0 0
## yellowblue 0 0 1587
## yellowgreen 0 0 0
## yelloworange 0 0 0
## yellowpink 0 0 0
## yellowpurple 0 0 0
## yellowred 0 0 0
# hist(unique(sl$IT_imputed[sl$trt2 == "low"]))
# rug(unique(sl$IT_imputed[sl$trt2 == "low"]))
#
# hist(unique(sl$IT_imputed[sl$trt2 == "high"]))
# rug(unique(sl$IT_imputed[sl$trt2 == "high"]))
# number of trials per bee and treatment
xtabs(~sl$beeCol + sl$trialNum)
## sl$trialNum
## sl$beeCol 1 2 3 4 5 6 7 8 9 10 11 12
## blue 36 0 0 0 0 0 0 0 0 0 0 0
## gold 54 0 0 0 0 0 0 0 0 0 0 0
## goldred 3 0 0 0 0 0 0 0 0 0 0 0
## green 9 4 0 0 0 0 0 0 0 0 0 0
## lime 28 0 0 0 0 0 0 0 0 0 0 0
## limeblue 50 54 50 50 53 55 52 62 53 51 0 0
## limegold 54 0 0 0 0 0 0 0 0 0 0 0
## limegreen 50 0 0 0 0 0 0 0 0 0 0 0
## limeorange 33 51 0 0 0 0 0 0 0 0 0 0
## limepink 51 0 0 0 0 0 0 0 0 0 0 0
## limepurple 53 243 56 0 0 0 0 0 0 0 0 0
## limepurpleyellow 58 0 0 0 0 0 0 0 0 0 0 0
## limered 50 281 351 298 486 537 393 170 96 105 707 0
## limesilver 3 0 0 0 0 0 0 0 0 0 0 0
## limewhite 50 51 0 0 0 0 0 0 0 0 0 0
## limeyellow 52 0 0 0 0 0 0 0 0 0 0 0
## orange 85 0 11 88 14 85 137 81 0 0 0 0
## orangeblue 50 0 0 0 0 0 0 0 0 0 0 0
## orangegreen 50 0 0 0 0 0 0 0 0 0 0 0
## orangepink 50 0 0 0 0 0 0 0 0 0 0 0
## orangepurple 50 0 0 0 0 0 0 0 0 0 0 0
## purple 9 0 0 0 0 0 0 0 0 0 0 0
## redblue 49 60 67 131 72 50 60 55 33 49 47 0
## redgreen 51 53 51 53 51 60 54 54 52 51 0 0
## redpink 56 58 71 57 62 62 76 119 62 54 0 0
## redpurple 55 90 73 0 0 0 0 0 0 0 0 0
## silver 26 0 0 0 0 0 0 0 0 0 0 0
## white 48 67 46 40 57 25 0 0 0 0 0 0
## whiteblue 56 314 307 190 539 537 105 103 276 111 248 0
## whitegold 48 66 76 174 179 198 216 135 64 69 0 0
## whitegreen 39 0 0 0 0 0 0 0 0 0 0 0
## whiteorange 54 260 432 376 0 0 0 0 0 0 0 0
## whitepink 64 143 71 99 78 165 122 127 105 118 92 0
## whitepurple 55 0 0 0 0 0 0 0 0 0 0 0
## whitered 59 195 169 278 109 89 112 87 138 108 0 0
## whiteyellow 55 47 55 0 0 0 0 0 0 0 0 0
## yellowblue 54 128 92 114 149 116 251 198 167 118 200 0
## yellowgreen 51 50 56 50 58 52 56 51 54 54 0 0
## yelloworange 54 176 205 475 113 256 111 141 252 110 220 0
## yellowpink 58 142 103 140 235 187 278 111 366 266 372 51
## yellowpurple 54 346 360 330 153 0 0 0 0 0 0 0
## yellowred 57 56 58 55 51 57 55 56 50 52 0 0
# tabulate num bees in each trt
tdf <- data.frame(xtabs(~sl$beeCol + sl$trt))
tdf2 <- acast(tdf, sl.beeCol~sl.trt, value.var="Freq")
tdf3 <- as.data.frame(tdf2 > 0)
trtdf<-table(interaction(tdf3$full, tdf3$high, tdf3$low))
names(trtdf) <- c("FULL ONLY", "FULL-HIGH", "FULL-LOW", "NA")
trtdf
## FULL ONLY FULL-HIGH FULL-LOW NA
## 28 6 8 0
# get number of buzzes in each treatment
colSums(tdf2)
## full high low
## 5095 6505 12703
# get the number of bees for each treatment (every bee has "full" treatment)
apply(tdf2, MARGIN = 2, function(x) length(x[x>0]))
## full high low
## 42 6 8
head(sl)
## beeCol hive meanFreq IT_imputed index freq amp
## 1 blue 4 395.2778 3.61 1 400 0.06328
## 2 blue 4 395.2778 3.61 2 360 0.34299
## 3 blue 4 395.2778 3.61 3 430 0.44099
## 4 blue 4 395.2778 3.61 4 420 0.16841
## 5 blue 4 395.2778 3.61 5 420 0.15284
## 6 blue 4 395.2778 3.61 6 410 0.18302
## datetime rewNum rewTF lowRewAmp highrewAmp
## 1 2016_11_22__08_23_47_721 1 T 0 5
## 2 2016_11_22__08_23_49_677 2 T 0 5
## 3 2016_11_22__08_23_50_139 3 T 0 5
## 4 2016_11_22__08_24_10_938 4 T 0 5
## 5 2016_11_22__08_24_18_013 5 T 0 5
## 6 2016_11_22__08_24_18_464 6 T 0 5
## BeeNumCol
## 1 BeeBlue_22Nov2016_Hive4_initial
## 2 BeeBlue_22Nov2016_Hive4_initial
## 3 BeeBlue_22Nov2016_Hive4_initial
## 4 BeeBlue_22Nov2016_Hive4_initial
## 5 BeeBlue_22Nov2016_Hive4_initial
## 6 BeeBlue_22Nov2016_Hive4_initial
## accFile trialNum
## 1 2016_11_22__08_23_47_721_220_450_test.txt 1
## 2 2016_11_22__08_23_49_677_220_450_test.txt 1
## 3 2016_11_22__08_23_50_139_220_450_test.txt 1
## 4 2016_11_22__08_24_10_938_220_450_test.txt 1
## 5 2016_11_22__08_24_18_013_220_450_test.txt 1
## 6 2016_11_22__08_24_18_464_220_450_test.txt 1
## datetime_str lowFrq highFrq IT trt amp_acc
## 1 2016-11-22 08:23:47.721000 220 450 3.61 full 6.222222
## 2 2016-11-22 08:23:49.677000 220 450 3.61 full 33.725664
## 3 2016-11-22 08:23:50.139000 220 450 3.61 full 43.361849
## 4 2016-11-22 08:24:10.938000 220 450 3.61 full 16.559489
## 5 2016-11-22 08:24:18.013000 220 450 3.61 full 15.028515
## 6 2016-11-22 08:24:18.464000 220 450 3.61 full 17.996067
dim(sl) # should be 24303 rows
## [1] 24303 21
# hive 5 is most common
table(sl$hive, useNA = 'always')
##
## 3 4 5 <NA>
## 513 1047 22743 0
## make sure all bee colors are lowercase
sl$beeCol <- tolower(sl$beeCol)
# make sure there are values lower than 220 and higher than 450
# (the cutoff for buzzes used in the experiment)
hist(sl$freq, breaks = seq(215, 450, by = 10))
nrow(sl[sl$freq < 220 | sl$freq > 450,]) # should have 0 rows
## [1] 0
# look at treatments
xtabs(~sl$beeCol+ trt, data = sl )
## trt
## sl$beeCol full high low
## blue 36 0 0
## gold 54 0 0
## goldred 3 0 0
## green 13 0 0
## lime 28 0 0
## limeblue 530 0 0
## limegold 54 0 0
## limegreen 50 0 0
## limeorange 84 0 0
## limepink 51 0 0
## limepurple 109 0 243
## limepurpleyellow 58 0 0
## limered 50 0 3424
## limesilver 3 0 0
## limewhite 101 0 0
## limeyellow 52 0 0
## orange 85 416 0
## orangeblue 50 0 0
## orangegreen 50 0 0
## orangepink 50 0 0
## orangepurple 50 0 0
## purple 9 0 0
## redblue 673 0 0
## redgreen 530 0 0
## redpink 56 0 621
## redpurple 55 163 0
## silver 26 0 0
## white 283 0 0
## whiteblue 56 0 2730
## whitegold 48 0 1177
## whitegreen 39 0 0
## whiteorange 54 0 1068
## whitepink 135 1049 0
## whitepurple 55 0 0
## whitered 59 1285 0
## whiteyellow 157 0 0
## yellowblue 54 1533 0
## yellowgreen 532 0 0
## yelloworange 54 2059 0
## yellowpink 58 0 2251
## yellowpurple 54 0 1189
## yellowred 547 0 0
# find percentage reward by treatment
mean(grepl("[tT]", as.character(sl$rewTF))) # overall mean
## [1] 0.4155454
# percentage that were rewarded by treatment
tapply((grepl("[tT]", as.character(sl$rewTF))), INDEX= sl$trt, mean)
## full high low
## 1.0000000 0.3535742 0.2128631
# total number of trials for each treatment
tapply((grepl("[tT]", as.character(sl$rewTF))), INDEX= sl$trt, length)
## full high low
## 5095 6505 12703
# total number of rewards per treatment
tapply((grepl("[tT]", as.character(sl$rewTF))), INDEX= sl$trt, FUN = function(x) sum(x))
## full high low
## 5095 2300 2704
# total number of trials that were unrewarded per treatment
tapply((grepl("[tT]", as.character(sl$rewTF))), INDEX= sl$trt, FUN = function(x) sum(!x))
## full high low
## 0 4205 9999
# calculate mean freq per bee per trial
d1 <- as.data.frame((tapply(sl$freq, INDEX = list(sl$beeCol, sl$trt), mean)))
d2 <- as.data.frame((tapply(sl$freq, INDEX = list(sl$beeCol, sl$trt), length)))
d3 <- data.frame(d1, d2)
d3$beeCol = row.names(d3)
d4 <- melt(d3, id.vars=c("beeCol"), measure.vars=c("full", "high", "low" ), value.name = "buzzFreq")
d5 <- melt(d3, id.vars=c("beeCol"), measure.vars=c("full.1", "high.1", "low.1" ), value.name = "NumBuzzes")
d6 <- data.frame(d4, d5)
d7 <- d6[!is.na(d6$buzzFreq),]
d7
## beeCol variable buzzFreq beeCol.1 variable.1
## 1 blue full 395.2778 blue full.1
## 2 gold full 348.3333 gold full.1
## 3 goldred full 406.6667 goldred full.1
## 4 green full 308.4615 green full.1
## 5 lime full 374.6429 lime full.1
## 6 limeblue full 304.1509 limeblue full.1
## 7 limegold full 344.6296 limegold full.1
## 8 limegreen full 370.8000 limegreen full.1
## 9 limeorange full 321.4286 limeorange full.1
## 10 limepink full 308.4314 limepink full.1
## 11 limepurple full 340.6422 limepurple full.1
## 12 limepurpleyellow full 336.7241 limepurpleyellow full.1
## 13 limered full 354.6000 limered full.1
## 14 limesilver full 306.6667 limesilver full.1
## 15 limewhite full 309.7030 limewhite full.1
## 16 limeyellow full 339.4231 limeyellow full.1
## 17 orange full 362.4706 orange full.1
## 18 orangeblue full 301.4000 orangeblue full.1
## 19 orangegreen full 296.2000 orangegreen full.1
## 20 orangepink full 286.6000 orangepink full.1
## 21 orangepurple full 329.0000 orangepurple full.1
## 22 purple full 382.2222 purple full.1
## 23 redblue full 341.2927 redblue full.1
## 24 redgreen full 326.9057 redgreen full.1
## 25 redpink full 288.3929 redpink full.1
## 26 redpurple full 311.4545 redpurple full.1
## 27 silver full 305.3846 silver full.1
## 28 white full 319.1873 white full.1
## 29 whiteblue full 318.2143 whiteblue full.1
## 30 whitegold full 315.2083 whitegold full.1
## 31 whitegreen full 300.7692 whitegreen full.1
## 32 whiteorange full 324.4444 whiteorange full.1
## 33 whitepink full 333.0370 whitepink full.1
## 34 whitepurple full 328.1818 whitepurple full.1
## 35 whitered full 288.8136 whitered full.1
## 36 whiteyellow full 316.8153 whiteyellow full.1
## 37 yellowblue full 313.7037 yellowblue full.1
## 38 yellowgreen full 284.7368 yellowgreen full.1
## 39 yelloworange full 268.1481 yelloworange full.1
## 40 yellowpink full 323.6207 yellowpink full.1
## 41 yellowpurple full 354.4444 yellowpurple full.1
## 42 yellowred full 337.4589 yellowred full.1
## 59 orange high 381.5385 orange high.1
## 68 redpurple high 336.4417 redpurple high.1
## 75 whitepink high 326.2917 whitepink high.1
## 77 whitered high 321.5019 whitered high.1
## 79 yellowblue high 309.4455 yellowblue high.1
## 81 yelloworange high 311.5590 yelloworange high.1
## 95 limepurple low 344.2798 limepurple low.1
## 97 limered low 364.7138 limered low.1
## 109 redpink low 290.2576 redpink low.1
## 113 whiteblue low 348.4505 whiteblue low.1
## 114 whitegold low 329.6092 whitegold low.1
## 116 whiteorange low 361.1610 whiteorange low.1
## 124 yellowpink low 346.0640 yellowpink low.1
## 125 yellowpurple low 361.1102 yellowpurple low.1
## NumBuzzes
## 1 36
## 2 54
## 3 3
## 4 13
## 5 28
## 6 530
## 7 54
## 8 50
## 9 84
## 10 51
## 11 109
## 12 58
## 13 50
## 14 3
## 15 101
## 16 52
## 17 85
## 18 50
## 19 50
## 20 50
## 21 50
## 22 9
## 23 673
## 24 530
## 25 56
## 26 55
## 27 26
## 28 283
## 29 56
## 30 48
## 31 39
## 32 54
## 33 135
## 34 55
## 35 59
## 36 157
## 37 54
## 38 532
## 39 54
## 40 58
## 41 54
## 42 547
## 59 416
## 68 163
## 75 1049
## 77 1285
## 79 1533
## 81 2059
## 95 243
## 97 3424
## 109 621
## 113 2730
## 114 1177
## 116 1068
## 124 2251
## 125 1189
# this plot shows how individual bees changed
ggplot(d7, aes(x = variable, y = buzzFreq)) +
geom_boxplot(aes(fill = variable == "full"), alpha = 0.2) +
geom_point(aes(color = variable == "full", size = NumBuzzes)) +
geom_line(aes(group = beeCol), alpha = 0.2) +
scale_fill_viridis(name = "hihi", guide = FALSE, discrete = TRUE, option = "A", end = 0.6) +
scale_color_viridis(name = "hihi", guide = FALSE, discrete = TRUE, option = "A", end = 0.6) +
labs(x = "treatment")
Visualizations and modeling
# calculate trial averages and plot
sl$colNum = paste(sl$beeCol, sl$hive, sep = "_")
sl$trt <- as.character(sl$trt)
sl$trt[sl$trt == "full" & sl$trialNum >1 ] <- "full_2"
aggdata <- aggregate(sl$freq, by=list(colNum = sl$colNum,trialNum =sl$trialNum, trt = sl$trt), FUN=mean, na.rm=TRUE)
colnames(aggdata)[colnames(aggdata) == "x"] = "freq"
aggdata_sd <- aggregate(sl$freq, by=list(colNum = sl$colNum,trialNum =sl$trialNum, trt = sl$trt), FUN=sd, na.rm=TRUE)
colnames(aggdata_sd)[colnames(aggdata_sd) == "x"] = "freq_sd"
aggdata = merge(aggdata, aggdata_sd)
#aggdata$trt = sapply(1:nrow(aggdata), FUN = function(x){sl[sl$colNum == aggdata$colNum[x] & sl$trialNum == aggdata$trialNum[x], "trt"][1]})
aggdata = aggdata[order(aggdata$colNum, aggdata$trialNum, decreasing = FALSE), ]
agg_sm = aggdata[aggdata$trialNum <= 2, ]
rownames(agg_sm) = 1:nrow(agg_sm)
agg_sm
## colNum trialNum trt freq freq_sd
## 1 blue_4 1 full 395.2778 39.38717
## 2 gold_3 1 full 348.3333 37.15089
## 3 goldred_4 1 full 406.6667 32.14550
## 4 green_4 1 full 318.8889 57.32461
## 5 green_4 2 full_2 285.0000 59.16080
## 6 lime_5 1 full 374.6429 36.36055
## 7 limeblue_5 1 full 302.4000 34.73280
## 8 limeblue_5 2 full_2 299.8148 39.02059
## 9 limegold_5 1 full 344.6296 44.07168
## 10 limegreen_5 1 full 370.8000 29.40498
## 11 limeorange_5 1 full 285.7576 35.26953
## 12 limeorange_5 2 full_2 344.5098 26.70683
## 13 limepink_5 1 full 308.4314 43.23760
## 14 limepurple_5 1 full 352.2642 38.66232
## 15 limepurple_5 2 low 344.2798 32.55948
## 16 limepurpleyellow_5 1 full 336.7241 50.76219
## 17 limered_5 1 full 354.6000 37.91559
## 18 limered_5 2 low 365.3381 36.20344
## 19 limesilver_5 1 full 306.6667 11.54701
## 20 limewhite_5 1 full 292.0000 44.90352
## 21 limewhite_5 2 full_2 327.0588 35.51305
## 22 limeyellow_5 1 full 339.4231 39.37818
## 23 orange_3 1 full 388.5294 38.22837
## 24 orange_5 1 full 345.0980 33.60789
## 25 orangeblue_5 1 full 301.4000 45.17585
## 26 orangegreen_5 1 full 296.2000 40.45103
## 27 orangepink_5 1 full 286.6000 27.15113
## 28 orangepurple_5 1 full 329.0000 25.65469
## 29 purple_3 1 full 382.2222 21.08185
## 30 redblue_4 1 full 342.2449 40.53113
## 31 redblue_4 2 full_2 335.6667 51.89189
## 32 redgreen_5 1 full 334.7059 34.19666
## 33 redgreen_5 2 full_2 314.5283 31.53545
## 34 redpink_5 1 full 288.3929 41.06654
## 35 redpink_5 2 low 275.3448 34.24450
## 36 redpurple_5 1 full 311.4545 31.64726
## 37 redpurple_5 2 high 336.4444 40.75995
## 38 silver_5 1 full 305.3846 38.39070
## 39 white_4 1 full 343.1250 47.49720
## 40 white_4 2 full_2 342.2388 38.95689
## 41 whiteblue_5 1 full 318.2143 27.17667
## 42 whiteblue_5 2 low 351.6242 29.13402
## 43 whitegold_5 1 full 315.2083 34.20772
## 44 whitegold_5 2 low 307.2727 45.76009
## 45 whitegreen_4 1 full 300.7692 33.43372
## 46 whiteorange_5 1 full 324.4444 51.71280
## 47 whiteorange_5 2 low 355.2308 44.74234
## 48 whitepink_5 1 full 336.8750 42.68285
## 49 whitepink_5 2 high 323.9161 28.55770
## 50 whitepurple_5 1 full 328.1818 43.50781
## 51 whitered_5 1 full 288.8136 36.05956
## 52 whitered_5 2 high 304.6667 42.10035
## 53 whiteyellow_5 1 full 303.4545 39.07344
## 54 whiteyellow_5 2 full_2 312.5532 36.56252
## 55 yellowblue_5 1 full 313.7037 44.86021
## 56 yellowblue_5 2 high 327.9688 49.76099
## 57 yellowgreen_5 1 full 298.0392 40.64577
## 58 yellowgreen_5 2 full_2 300.6000 27.58364
## 59 yelloworange_5 1 full 268.1481 36.18996
## 60 yelloworange_5 2 high 313.9205 39.62646
## 61 yellowpink_5 1 full 323.6207 39.98676
## 62 yellowpink_5 2 low 335.7042 41.00507
## 63 yellowpurple_5 1 full 354.4444 36.37678
## 64 yellowpurple_5 2 low 373.0058 34.17237
## 65 yellowred_5 1 full 320.1754 41.85396
## 66 yellowred_5 2 full_2 333.0357 29.96481
agg_sm[duplicated(data.frame(agg_sm$colNum, agg_sm$trialNum)), ]
## [1] colNum trialNum trt freq freq_sd
## <0 rows> (or 0-length row.names)
diffdf <- sapply(unique(agg_sm$colNum), FUN = function(x){
tmp = agg_sm[agg_sm$colNum == x, ]
if(nrow(tmp) <= 1)
diff = NA
else
diff = tmp$freq[tmp$trialNum == 2] - tmp$freq[tmp$trialNum == 1]
return(diff)
})
trtDF = sapply(unique(agg_sm$colNum), FUN = function(x){
tmp = agg_sm[agg_sm$colNum == x, ]
ttrs = paste(tmp$trt[tmp$trialNum == 1], tmp$trt[tmp$trialNum == 2], sep = "_")
return(ttrs)
})
buzzdiffs = data.frame(trtDF, diffdf)
ggplot(buzzdiffs, aes(x = trtDF, y= diffdf)) +
geom_boxplot() +
geom_point()
## Warning: Removed 20 rows containing non-finite values (stat_boxplot).
## Warning: Removed 20 rows containing missing values (geom_point).
agg2 = aggregate(sl$freq, by=list(colNum = sl$colNum, fullTrt = sl$trt == "full", trt = sl$trt), FUN=mean, na.rm=TRUE)
colnames(agg2)[colnames(agg2) == "x"] = "freq"
agg2$trt = as.character(agg2$trt)
diffdf <- t(as.data.frame(t(sapply(unique(agg2$colNum), FUN = function(x){
tmp = agg2[agg2$colNum == x, ]
if(nrow(tmp) <= 1)
return(NA)
if (length(unique(tmp$trt)) > 2){
tmp = tmp[tmp$trt != "full_2", ]
}
diff = tmp$freq[!tmp$fullTrt] - tmp$freq[tmp$fullTrt]
return(diff)
}))))
length(diffdf)
## [1] 43
trtDF = sapply(unique(agg2$colNum), FUN = function(x){
tmp = agg2[agg2$colNum == x, ]
if (length(unique(tmp$trt)) > 2){
tmp = tmp[tmp$trt != "full_2", ]
}
ttrs = paste(tmp$trt[tmp$fullTrt], tmp$trt[!tmp$fullTrt], sep = "_")
return(ttrs)
})
length(trtDF)
## [1] 43
buzzdiffs = data.frame(trtDF, diffdf)
tapply(buzzdiffs$diffdf, INDEX = buzzdiffs$trtDF, mean)
## full_ full_full_2 full_high full_low
## NA 4.851689 13.209001 14.307123
ggplot(droplevels(buzzdiffs[buzzdiffs$trtDF != "full_", ]), aes(x = trtDF, y= as.numeric(diffdf))) +
geom_boxplot() +
geom_point() +
labs(y = "diff in frequency -- averages")
MODELING
sl$trt = relevel(factor(sl$trt), ref = "full")
sl$trt2 = mapvalues(sl$trt, from = c("full", "full_2", "high", "low"),
to = c("full", "full", "high", "low"))
# summary for paper
# fit a varying slope and intercept for colNum (bee ID), and allow the slope of the trialNum
# variable to vary by colNum (beeID)
m2 = lmer(freq ~ trt2 + hive + trialNum + IT_imputed + (1+trialNum|colNum), data = sl, REML = FALSE)
# this model estimates a global intercept
# random effect intercept for colNum (beeID)
# a single global estimate for trialNum
# the effect of trialNum within each level of colNum
# the correlation between intercept of trialNum across levels of colNum
vif.mer <- function (fit) {
## adapted from rms::vif
v <- vcov(fit)
nam <- names(fixef(fit))
## exclude intercepts
ns <- sum(1 * (nam == "Intercept" | nam == "(Intercept)"))
if (ns > 0) {
v <- v[-(1:ns), -(1:ns), drop = FALSE]
nam <- nam[-(1:ns)]
}
d <- diag(v)^0.5
v <- diag(solve(v/(d %o% d)))
names(v) <- nam
v
}
vif.mer(m2) # fine
## trt2high trt2low hive4 hive5 trialNum IT_imputed
## 1.010841 1.004759 2.648702 2.664859 1.014513 1.143183
m2_1 = lmer(freq ~ trt2* IT_imputed + hive + trialNum + (1+trialNum|colNum), data = sl, REML = FALSE)
# compare models with BIC
BIC(m2, m2_1) # m2 (no interaction) is better
## df BIC
## m2 11 246086.5
## m2_1 13 246090.2
anova(m2, m2_1) # Note: this disagrees with BIC
## Data: sl
## Models:
## m2: freq ~ trt2 + hive + trialNum + IT_imputed + (1 + trialNum |
## m2: colNum)
## m2_1: freq ~ trt2 * IT_imputed + hive + trialNum + (1 + trialNum |
## m2_1: colNum)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2 11 245997 246087 -122988 245975
## m2_1 13 245985 246090 -122979 245959 16.547 2 0.0002552 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m3 = update(m2, .~. - hive)
BIC(m2, m3) # m3 is better (without hive)
## df BIC
## m2 11 246086.5
## m3 9 246073.9
anova(m2, m3) # Note: disagrees with BIC
## Data: sl
## Models:
## m3: freq ~ trt2 + trialNum + IT_imputed + (1 + trialNum | colNum)
## m2: freq ~ trt2 + hive + trialNum + IT_imputed + (1 + trialNum |
## m2: colNum)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m3 9 246001 246074 -122991 245983
## m2 11 245997 246087 -122988 245975 7.5417 2 0.02303 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# refit best model with REML = TRUE
m3 <- update(m3, .~., REML = TRUE)
summary(m3) # summary for paper -- leave trial number in, to account for random effect
## Linear mixed model fit by REML ['lmerMod']
## Formula: freq ~ trt2 + trialNum + IT_imputed + (1 + trialNum | colNum)
## Data: sl
##
## REML criterion at convergence: 245962
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0242 -0.5161 0.1663 0.6809 3.9135
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## colNum (Intercept) 1072 32.74
## trialNum 112 10.58 -0.73
## Residual 1439 37.93
## Number of obs: 24303, groups: colNum, 43
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 479.138 44.361 10.801
## trt2high 11.161 2.294 4.865
## trt2low 17.718 1.916 9.248
## trialNum 1.448 2.133 0.679
## IT_imputed -37.291 10.773 -3.462
##
## Correlation of Fixed Effects:
## (Intr) trt2hg trt2lw trilNm
## trt2high 0.069
## trt2low -0.005 0.009
## trialNum 0.007 -0.042 -0.037
## IT_imputed -0.993 -0.072 0.000 -0.093
plot(m3)
qqnorm(ranef(m3)$colNum[[1]])
qqline(ranef(m3)$colNum[[1]])
sjp.lmer(m3) # plot random effects to find any outliers
## Plotting random effects...
sjp.lmer(m3, type = 'fe', sort = TRUE, p.kr = FALSE) # plot fixed effects
## Computing p-values via Wald-statistics approximation (treating t as Wald z).
# post-hoc tests -- pvals for paper
summary(glht(m3, linfct = mcp(trt2 = "Tukey")), test = adjusted("none"))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = freq ~ trt2 + trialNum + IT_imputed + (1 + trialNum |
## colNum), data = sl, REML = TRUE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## high - full == 0 11.161 2.294 4.865 1.14e-06 ***
## low - full == 0 17.718 1.916 9.248 < 2e-16 ***
## low - high == 0 6.557 2.975 2.204 0.0275 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- none method)
## bonf adjusted pvals
summary(glht(m3, linfct = mcp(trt2 = "Tukey")), test = adjusted("bonf"))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = freq ~ trt2 + trialNum + IT_imputed + (1 + trialNum |
## colNum), data = sl, REML = TRUE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## high - full == 0 11.161 2.294 4.865 3.43e-06 ***
## low - full == 0 17.718 1.916 9.248 < 2e-16 ***
## low - high == 0 6.557 2.975 2.204 0.0826 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)
Generate CI’s
# set number of bOotstrap samples
nbootSims = 1000
ITmean = mean(tapply(sl$IT_imputed, INDEX = sl$beeCol, FUN = function(x) x[1] ))
print(ITmean)
## [1] 4.092619
# don't need hive, because that's not in the model we chose (above)
pframe <- data.frame(trt2 = levels(droplevels(sl$trt2)), IT_imputed = ITmean, colNum = 99999, trialNum = 2)
pframe$freq <- 0
pp <- predict(m3, newdata = pframe, re.form=NA, type = 'response') # re.form sets all random effects to 0
### Calculate CI's (using bootstrap, not accounting for random effects)
bb2 <- bootMer(m3, FUN=function(x) predict(x, pframe, re.form=NA, type = 'response'), nsim = nbootSims)
print(paste("Number of bootstrap samples", nrow(bb2$t)))
## [1] "Number of bootstrap samples 1000"
bb2_se <-apply(bb2$t,2,function(x) quantile(x, probs = c(0.025, 0.975)))
pframe$blo<-bb2_se[1,]
pframe$bhi<-bb2_se[2,]
pframe$predMean <- pp
pframe <- pframe[, c('trt2', "blo", "bhi", "predMean")]
pframe
## trt2 blo bhi predMean
## 1 full 321.7666 336.2141 329.4150
## 2 high 332.6790 348.3194 340.5761
## 3 low 339.1716 355.1462 347.1330
Make plots for paper
#plot 95% confidence intervals
# "Mean and bootstrap CI based on fixed-effects uncertainty ONLY"
# rename labels for plot
pframe$trt3 = mapvalues(pframe$trt2, from = c("full", "high", "low"),
to = c("Full range\n(220 - 450 Hz)", "High range\n(340 - 390 Hz)", "Low range\n(220 - 330 Hz)"))
pframe
## trt2 blo bhi predMean trt3
## 1 full 321.7666 336.2141 329.4150 Full range\n(220 - 450 Hz)
## 2 high 332.6790 348.3194 340.5761 High range\n(340 - 390 Hz)
## 3 low 339.1716 355.1462 347.1330 Low range\n(220 - 330 Hz)
g0 <- ggplot(pframe, aes(x=trt3, y=predMean))+
geom_point()+
labs(y = "Sonication frequency (Hz)", x = "Frequency range for reward") +
geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1)+
theme(plot.background = element_rect(fill = "transparent",colour = NA),
panel.background = element_rect(fill = "transparent",colour = NA)) +
annotate(geom="text", x=c(1,2,3), y=c(0, 0, 0) + 360, label=c("a", "b", "c"),
color="black")
g0
ggsave(plot = g0, filename = file.path(figDir, "SonicationFreqPredsAndCI_unadjusted.png"), width = 5, height = 4)
# save with transparent background
g0 <- ggplot(pframe, aes(x=trt3, y=predMean))+
geom_point()+
labs(y = "Sonication frequency (Hz)", x = "Frequency range for reward") +
geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1)+
theme(plot.background = element_rect(fill = "transparent",colour = NA),
panel.background = element_rect(fill = "transparent",colour = NA)) +
annotate(geom="text", x=c(1,2,3), y=c(0, 0, 0) + 360, label=c("a", "b", "b"),
color="black")
g0
ggsave(plot = g0, filename = file.path(figDir, "SonicationFreqPredsAndCI_adjustedPvals.pdf"), width = 5, height = 4)
# use the "effects" library to get CI's
# this basically gives the same result as bootstrapping
# it's a good double check
ee <- as.data.frame(Effect(c("trt2"),
fixed.predictors =list(given.values = c(trialNum = 2,
IT_imputed = ITmean)),
m3) )
# compare two methods -- very similar
ee
## trt2 fit se lower upper
## 1 full 329.4150 3.798940 321.9689 336.8612
## 2 high 340.5761 4.280727 332.1856 348.9666
## 3 low 347.1330 4.081495 339.1330 355.1330
pframe
## trt2 blo bhi predMean trt3
## 1 full 321.7666 336.2141 329.4150 Full range\n(220 - 450 Hz)
## 2 high 332.6790 348.3194 340.5761 High range\n(340 - 390 Hz)
## 3 low 339.1716 355.1462 347.1330 Low range\n(220 - 330 Hz)
g1 <- ggplot(ee, aes(x=trt2, y=fit))+
geom_point()+
labs(y = "Sonication Frequency (Hz)", x = "Treatment") +
geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.1)+
theme_classic() +
annotate(geom="text", x=c(1,2,3), y=c(0, 0, 0) + 360, label=c("a", "b", "c"),
color="black")
g1
____ END OF ANALYSIS FOR PAPER ____
Analysis for amplitude
# refref: here try log amp_acc and REML = FALSE with BIC
# interaction model
# note log-transformation to make model fit assumptions better
maa0 = lmer(log(amp_acc) ~ trt2* IT_imputed + hive + trialNum + (1+trialNum|colNum), data = sl, REML = FALSE)
# main effect model
maa1 = lmer(log(amp_acc) ~ trt2 + hive + trialNum + IT_imputed + (1+trialNum|colNum), data = sl, REML = FALSE)
BIC(maa0, maa1) # use no interaction (keep maa1)
## df BIC
## maa0 13 44377.70
## maa1 11 44360.43
anova(maa0, maa1) # note that this agrees with likelihood ratio test
## Data: sl
## Models:
## maa1: log(amp_acc) ~ trt2 + hive + trialNum + IT_imputed + (1 + trialNum |
## maa1: colNum)
## maa0: log(amp_acc) ~ trt2 * IT_imputed + hive + trialNum + (1 + trialNum |
## maa0: colNum)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## maa1 11 44271 44360 -22125 44249
## maa0 13 44272 44378 -22123 44246 2.9209 2 0.2321
maa2 = update(maa1, .~. - trt2)
BIC(maa1, maa2) # keep treatment (maa1)
## df BIC
## maa1 11 44360.43
## maa2 9 44375.90
anova(maa1, maa2) # agrees with LRT; p-val for paper, if needed
## Data: sl
## Models:
## maa2: log(amp_acc) ~ hive + trialNum + IT_imputed + (1 + trialNum |
## maa2: colNum)
## maa1: log(amp_acc) ~ trt2 + hive + trialNum + IT_imputed + (1 + trialNum |
## maa1: colNum)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## maa2 9 44303 44376 -22142 44285
## maa1 11 44271 44360 -22125 44249 35.665 2 1.801e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
maa3 <- update(maa1, .~. - hive)
BIC(maa1, maa3) # remove hive (maa3)
## df BIC
## maa1 11 44360.43
## maa3 9 44345.85
anova(maa1, maa3) # agrees with LRT
## Data: sl
## Models:
## maa3: log(amp_acc) ~ trt2 + trialNum + IT_imputed + (1 + trialNum |
## maa3: colNum)
## maa1: log(amp_acc) ~ trt2 + hive + trialNum + IT_imputed + (1 + trialNum |
## maa1: colNum)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## maa3 9 44273 44346 -22128 44255
## maa1 11 44271 44360 -22125 44249 5.6143 2 0.06038 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# refit with REML = TRUE for paper
maa3 <- update(maa3, .~., REML = TRUE)
summary(maa3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log(amp_acc) ~ trt2 + trialNum + IT_imputed + (1 + trialNum |
## colNum)
## Data: sl
##
## REML criterion at convergence: 44279.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.8377 -0.5961 0.1156 0.7367 2.9119
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## colNum (Intercept) 0.07079 0.26607
## trialNum 0.00203 0.04506 -0.57
## Residual 0.35872 0.59893
## Number of obs: 24303, groups: colNum, 43
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.398060 0.455225 5.268
## trt2high 0.194526 0.034446 5.647
## trt2low 0.067759 0.029326 2.311
## trialNum 0.005513 0.010107 0.545
## IT_imputed 0.324122 0.110268 2.939
##
## Correlation of Fixed Effects:
## (Intr) trt2hg trt2lw trilNm
## trt2high 0.107
## trt2low -0.009 0.029
## trialNum 0.036 -0.089 -0.083
## IT_imputed -0.995 -0.114 0.000 -0.087
# diagnostics
plot(maa3)
qqnorm(ranef(maa1)$colNum[[1]])
qqline(ranef(maa1)$colNum[[1]])
sjp.lmer(maa3) # plot random effects to find any outliers
## Plotting random effects...
sjp.lmer(maa3, type = 'fe', sort = TRUE, p.kr = FALSE) # plot fixed effects
## Computing p-values via Wald-statistics approximation (treating t as Wald z).
# post-hoc tests -- pvals for paper
summary(glht(maa3, linfct = mcp(trt2 = "Tukey")), test = adjusted("none"))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = log(amp_acc) ~ trt2 + trialNum + IT_imputed +
## (1 + trialNum | colNum), data = sl, REML = TRUE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## high - full == 0 0.19453 0.03445 5.647 1.63e-08 ***
## low - full == 0 0.06776 0.02933 2.311 0.02086 *
## low - high == 0 -0.12677 0.04459 -2.843 0.00447 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- none method)
# post-hoc tests with bonf adjustment
summary(glht(maa3, linfct = mcp(trt2 = "Tukey")), test = adjusted("bonf"))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = log(amp_acc) ~ trt2 + trialNum + IT_imputed +
## (1 + trialNum | colNum), data = sl, REML = TRUE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## high - full == 0 0.19453 0.03445 5.647 4.89e-08 ***
## low - full == 0 0.06776 0.02933 2.311 0.0626 .
## low - high == 0 -0.12677 0.04459 -2.843 0.0134 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)
Generate CI’s for amplitude
# set number of bootstrap samples
nbootSims2 = 1000
# don't need to include hive, but will set trialNum to 2
pframe <- data.frame(trt2 = levels(droplevels(sl$trt2)),
IT_imputed = ITmean, colNum = 99999, trialNum = 2)
pframe$amp_acc <- 0
# exponentiate to get back to original scale
pp <- exp(predict(maa3, newdata = pframe, re.form=NA, type = 'response')) # re.form sets all random effects to 0
### Calculate CI's (using bootstrap, not accounting for random effects)
bb2 <- bootMer(maa3, FUN=function(x) predict(x, pframe, re.form=NA, type = 'response'), nsim = nbootSims2)
print(paste("Number of bootstrap samples", nrow(bb2$t)))
## [1] "Number of bootstrap samples 1000"
bb2_se <-apply(bb2$t,2,function(x) quantile(x, probs = c(0.025, 0.975)))
#exponentiate to get back to orignal scale
pframe$blo<-exp(bb2_se[1,])
pframe$bhi<-exp(bb2_se[2,])
pframe$predMean <- pp
pframe <- pframe[, c('trt2', "blo", "bhi", "predMean")]
pframe
## trt2 blo bhi predMean
## 1 full 39.05177 45.01982 41.91285
## 2 high 46.44376 56.04234 50.91299
## 3 low 41.17250 48.98075 44.85125
Make plots for amplitude for paper
#plot 95% confidence intervals
# "Mean and bootstrap CI based on fixed-effects uncertainty ONLY"
pframe$trt3 = mapvalues(pframe$trt2, from = c("full", "high", "low"),
to = c("Full range\n(220 - 450 Hz)", "High range\n(340 - 390 Hz)", "Low range\n(220 - 330 Hz)"))
pframe
## trt2 blo bhi predMean trt3
## 1 full 39.05177 45.01982 41.91285 Full range\n(220 - 450 Hz)
## 2 high 46.44376 56.04234 50.91299 High range\n(340 - 390 Hz)
## 3 low 41.17250 48.98075 44.85125 Low range\n(220 - 330 Hz)
ga0 <- ggplot(pframe, aes(x=trt3, y=predMean))+
geom_point()+
labs( y = expression ("Sonication amplitude "(m~s^{-2})), x = "Frequency range for reward") +
geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1)+
theme(axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = 'none') +
theme_classic() +
annotate(geom="text", x=c(1,2,3), y=c(0, 0, 0) + 58, label=c("a", "b", "c"),
color="black")
ga0
ggsave(plot = ga0, filename = file.path(figDir, "SonicationAmpPredsAndCI_unadjusted.pdf"), width = 5, height = 4)
# make plot using bonferroni adjustments
ga0 <- ggplot(pframe, aes(x=trt3, y=predMean))+
geom_point()+
labs( y = expression ("Sonication amplitude "(m~s^{-2})), x = "Frequency range for reward") +
geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1)+
theme(axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = 'none') +
theme_classic() +
annotate(geom="text", x=c(1,2,3), y=c(0, 0, 0) + 58, label=c("a", "b", "a"),
color="black")
ga0
ggsave(plot = ga0, filename = file.path(figDir, "SonicationAmpPredsAndCI_BonfAdjusted.pdf"), width = 5, height = 4)
Visualize aggregated data for amplitude
aggdata <- aggregate(sl$amp_acc, by=list(colNum = sl$colNum,trialNum =sl$trialNum, trt = sl$trt), FUN=mean, na.rm=TRUE)
colnames(aggdata)[colnames(aggdata) == "x"] = "amp"
aggdata_sd <- aggregate(sl$amp_acc, by=list(colNum = sl$colNum,trialNum =sl$trialNum, trt = sl$trt), FUN=sd, na.rm=TRUE)
colnames(aggdata_sd)[colnames(aggdata_sd) == "x"] = "amp_sd"
aggdata = merge(aggdata, aggdata_sd)
aggdata = aggdata[order(aggdata$colNum, aggdata$trialNum, decreasing = FALSE), ]
agg_sm = aggdata[aggdata$trialNum <= 2, ]
rownames(agg_sm) = 1:nrow(agg_sm)
agg_sm
## colNum trialNum trt amp amp_sd
## 1 blue_4 1 full 40.03641 27.005004
## 2 gold_3 1 full 32.58012 13.395404
## 3 goldred_4 1 full 27.49328 8.575236
## 4 green_4 1 full 44.14520 45.654277
## 5 green_4 2 full_2 19.19985 21.543257
## 6 lime_5 1 full 72.53866 33.940021
## 7 limeblue_5 1 full 49.20279 24.837998
## 8 limeblue_5 2 full_2 59.04592 36.444346
## 9 limegold_5 1 full 34.95839 26.390545
## 10 limegreen_5 1 full 52.56061 32.949195
## 11 limeorange_5 1 full 43.28056 22.966005
## 12 limeorange_5 2 full_2 50.75113 25.907474
## 13 limepink_5 1 full 49.57493 19.711454
## 14 limepurple_5 1 full 43.98249 22.702402
## 15 limepurple_5 2 low 39.95722 19.359774
## 16 limepurpleyellow_5 1 full 34.54496 17.151702
## 17 limered_5 1 full 42.27701 20.963674
## 18 limered_5 2 low 40.63711 20.895896
## 19 limesilver_5 1 full 26.24910 7.336471
## 20 limewhite_5 1 full 37.65233 22.862123
## 21 limewhite_5 2 full_2 40.91118 19.936507
## 22 limeyellow_5 1 full 43.51579 23.156058
## 23 orange_3 1 full 32.25152 18.932874
## 24 orange_5 1 full 55.08235 20.197053
## 25 orangeblue_5 1 full 51.97426 23.301810
## 26 orangegreen_5 1 full 58.59866 33.410963
## 27 orangepink_5 1 full 50.98635 28.370137
## 28 orangepurple_5 1 full 34.53312 16.739467
## 29 purple_3 1 full 53.61466 21.935890
## 30 redblue_4 1 full 28.88506 11.312274
## 31 redblue_4 2 full_2 48.61796 25.277613
## 32 redgreen_5 1 full 78.72927 46.273733
## 33 redgreen_5 2 full_2 70.38830 47.080737
## 34 redpink_5 1 full 74.20294 37.536798
## 35 redpink_5 2 low 56.85081 31.951114
## 36 redpurple_5 1 full 58.06343 40.818616
## 37 redpurple_5 2 high 53.13637 26.813039
## 38 silver_5 1 full 46.34517 30.280907
## 39 white_4 1 full 47.76100 21.162920
## 40 white_4 2 full_2 46.45015 22.564532
## 41 whiteblue_5 1 full 66.11499 33.532056
## 42 whiteblue_5 2 low 67.00735 29.330027
## 43 whitegold_5 1 full 43.10847 26.148044
## 44 whitegold_5 2 low 45.20369 27.241996
## 45 whitegreen_4 1 full 38.87762 18.184245
## 46 whiteorange_5 1 full 61.63802 39.835205
## 47 whiteorange_5 2 low 45.85399 25.838373
## 48 whitepink_5 1 full 63.53760 35.030181
## 49 whitepink_5 2 high 80.19250 47.237615
## 50 whitepurple_5 1 full 47.77247 35.753931
## 51 whitered_5 1 full 48.37233 25.860700
## 52 whitered_5 2 high 74.79446 45.764818
## 53 whiteyellow_5 1 full 65.50425 32.808944
## 54 whiteyellow_5 2 full_2 56.00923 33.372943
## 55 yellowblue_5 1 full 52.45834 30.688932
## 56 yellowblue_5 2 high 79.36156 41.580725
## 57 yellowgreen_5 1 full 61.58710 33.835105
## 58 yellowgreen_5 2 full_2 79.18968 27.697874
## 59 yelloworange_5 1 full 73.32616 36.109449
## 60 yelloworange_5 2 high 54.48161 24.329362
## 61 yellowpink_5 1 full 52.74199 24.869065
## 62 yellowpink_5 2 low 57.76271 29.939244
## 63 yellowpurple_5 1 full 62.77375 29.437162
## 64 yellowpurple_5 2 low 68.98565 35.865129
## 65 yellowred_5 1 full 62.77866 39.753892
## 66 yellowred_5 2 full_2 47.97187 28.376661
agg_sm[duplicated(data.frame(agg_sm$colNum, agg_sm$trialNum)), ]
## [1] colNum trialNum trt amp amp_sd
## <0 rows> (or 0-length row.names)
ggplot(agg_sm, aes(x = trt, y = amp, fill = trialNum > 1)) +
geom_boxplot(alpha = 0.2) +
geom_point(aes(color = trialNum>1)) +
geom_line(aes(group = colNum))
diffdf <- sapply(unique(agg_sm$colNum), FUN = function(x){
tmp = agg_sm[agg_sm$colNum == x, ]
if(nrow(tmp) <= 1)
diff = NA
else
diff = tmp$amp[tmp$trialNum == 2] - tmp$amp[tmp$trialNum == 1]
return(diff)
})
trtDF = sapply(unique(agg_sm$colNum), FUN = function(x){
tmp = agg_sm[agg_sm$colNum == x, ]
ttrs = paste(tmp$trt[tmp$trialNum == 1], tmp$trt[tmp$trialNum == 2], sep = "_")
return(ttrs)
})
buzzdiffs = data.frame(trtDF, diffdf)
ggplot(buzzdiffs, aes(x = trtDF, y= diffdf)) +
geom_boxplot() +
geom_point()+
labs(y = "amplitude difference m/s/s")
## Warning: Removed 20 rows containing non-finite values (stat_boxplot).
## Warning: Removed 20 rows containing missing values (geom_point).
agg2 = aggregate(sl$amp_acc, by=list(colNum = sl$colNum, fullTrt = sl$trt == "full", trt = sl$trt), FUN=mean, na.rm=TRUE)
colnames(agg2)[colnames(agg2) == "x"] = "amp"
agg2$trt = as.character(agg2$trt)
diffdf <- t(as.data.frame(t(sapply(unique(agg2$colNum), FUN = function(x){
tmp = agg2[agg2$colNum == x, ]
if(nrow(tmp) <= 1)
return(NA)
if (length(unique(tmp$trt)) > 2){
tmp = tmp[tmp$trt != "full_2", ]
}
diff = tmp$amp[!tmp$fullTrt] - tmp$amp[tmp$fullTrt]
return(diff)
}))))
length(diffdf)
## [1] 43
trtDF = sapply(unique(agg2$colNum), FUN = function(x){
tmp = agg2[agg2$colNum == x, ]
if (length(unique(tmp$trt)) > 2){
tmp = tmp[tmp$trt != "full_2", ]
}
ttrs = paste(tmp$trt[tmp$fullTrt], tmp$trt[!tmp$fullTrt], sep = "_")
return(ttrs)
})
length(trtDF)
## [1] 43
buzzdiffs = data.frame(trtDF, diffdf)
tapply(buzzdiffs$diffdf, INDEX = buzzdiffs$trtDF, mean)
## full_ full_full_2 full_high full_low
## NA 2.871714 9.194234 3.249068
ggplot(droplevels(buzzdiffs[buzzdiffs$trtDF != "full_", ]), aes(x = trtDF, y= as.numeric(diffdf))) +
geom_boxplot() +
geom_point()+
labs(y = "amplitude difference m/s/s")
Plot amplitude vs. frequency for initial trial, accounting for size
ggplot(sl[sl$trialNum == 1, ], aes(x = freq, y = amp_acc)) +
geom_point(position = position_jitter(width = 3), alpha = 0.2) +
theme_classic() +
geom_smooth(aes(group = trt2)) +
labs( y = expression ("Sonication amplitude "(m~s^{-2})), x = "Frequency (Hz)")
## `geom_smooth()` using method = 'gam'
centerFreq = scale(sl[sl$trialNum == 1, "freq"])
m1 <- lmer(log(amp_acc)~centerFreq + I(centerFreq^2) + I(centerFreq^3) + IT_imputed + (1|beeCol), data = sl[sl$trialNum == 1, ])
summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log(amp_acc) ~ centerFreq + I(centerFreq^2) + I(centerFreq^3) +
## IT_imputed + (1 | beeCol)
## Data: sl[sl$trialNum == 1, ]
##
## REML criterion at convergence: 3526.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5361 -0.6344 0.0693 0.6857 2.4239
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeCol (Intercept) 0.08395 0.2897
## Residual 0.32877 0.5734
## Number of obs: 1971, groups: beeCol, 42
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.00589 0.57016 3.518
## centerFreq 0.37124 0.03004 12.357
## I(centerFreq^2) -0.08590 0.01227 -7.003
## I(centerFreq^3) -0.03835 0.01055 -3.636
## IT_imputed 0.43530 0.13776 3.160
##
## Correlation of Fixed Effects:
## (Intr) cntrFr I(F^2) I(F^3)
## centerFreq -0.047
## I(cntrFr^2) -0.093 0.081
## I(cntrFr^3) -0.018 -0.844 0.065
## IT_imputed -0.996 0.043 0.071 0.019
plot(m1)
# predict -- using mean IT span
ppdf <- data.frame(centerFreq = sort(unique(centerFreq)), IT_imputed = ITmean, beeCol = 99999, acc_amp = 0)
# exponentiate to get back to original scale
ppdf$predAmp = exp(predict(m1, newdata = ppdf, type = "response", re.form = NA))
cf_Unscaled = ppdf$centerFreq * attr(centerFreq, 'scaled:scale') + attr(centerFreq, 'scaled:center')
ggplot(ppdf, aes(x = cf_Unscaled, y = predAmp))+
geom_line() +
labs(x = "Sonication Frequency (Hz)", y = expression ("Sonication amplitude "(m~s^{-2})) ) +
geom_point(data = sl[sl$trialNum == 1, ],
aes(x = freq, y = amp_acc),
alpha = 0.2, position = position_jitter(width =2), pch =16, size = .5)
ggsave(filename = file.path(figDir, "FreqVsAmp_1stTrial_rawData.pdf"), width = 4, height = 3)
g22 <- ggplot(ppdf, aes(x = cf_Unscaled, y = predAmp))+
geom_line() +
labs(x = "Sonication Frequency (Hz)", y = expression ("Predicted sonication amplitude "(m~s^{-2})) )
g22
ggsave(filename = file.path(figDir, "FreqVsAmp_1stTrial.pdf"), width = 4, height = 3)
sessionInfo()
## R version 3.4.2 (2017-09-28)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.2
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] viridis_0.4.0 viridisLite_0.2.0 effects_4.0-0
## [4] carData_3.0-0 plyr_1.8.4 multcomp_1.4-8
## [7] TH.data_1.0-8 MASS_7.3-47 survival_2.41-3
## [10] mvtnorm_1.0-6 sjPlot_2.4.0 lme4_1.1-14
## [13] Matrix_1.2-11 reshape2_1.4.2 ggplot2_2.2.1
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-131 pbkrtest_0.4-7 RColorBrewer_1.1-2
## [4] rprojroot_1.2 tools_3.4.2 TMB_1.7.11
## [7] backports_1.1.1 R6_2.2.2 sjlabelled_1.0.5
## [10] DT_0.2 mgcv_1.8-20 lazyeval_0.2.1
## [13] colorspace_1.3-2 nnet_7.3-12 tidyselect_0.2.3
## [16] gridExtra_2.3 mnormt_1.5-5 compiler_3.4.2
## [19] sandwich_2.4-0 labeling_0.3 scales_0.5.0
## [22] lmtest_0.9-35 psych_1.7.8 blme_1.0-4
## [25] stringr_1.2.0 digest_0.6.12 foreign_0.8-69
## [28] minqa_1.2.4 rmarkdown_1.8 stringdist_0.9.4.6
## [31] pkgconfig_2.0.1 htmltools_0.3.6 pwr_1.2-1
## [34] htmlwidgets_0.9 rlang_0.1.4 shiny_1.0.5
## [37] bindr_0.1 zoo_1.8-0 dplyr_0.7.4
## [40] magrittr_1.5 modeltools_0.2-21 bayesplot_1.4.0
## [43] Rcpp_0.12.13 munsell_0.4.3 abind_1.4-5
## [46] prediction_0.2.0 stringi_1.1.6 yaml_2.1.14
## [49] merTools_0.3.0 snakecase_0.5.1 grid_3.4.2
## [52] parallel_3.4.2 sjmisc_2.6.2 forcats_0.2.0
## [55] lattice_0.20-35 ggeffects_0.2.2 haven_1.1.0
## [58] splines_3.4.2 sjstats_0.12.0 knitr_1.17
## [61] codetools_0.2-15 stats4_3.4.2 glue_1.2.0
## [64] evaluate_0.10.1 modelr_0.1.1 nloptr_1.0.4
## [67] httpuv_1.3.5 gtable_0.2.0 purrr_0.2.4
## [70] tidyr_0.7.2 assertthat_0.2.0 mime_0.5
## [73] coin_1.2-1 xtable_1.8-2 broom_0.4.2
## [76] survey_3.32-1 coda_0.19-1 tibble_1.3.4
## [79] arm_1.9-3 glmmTMB_0.1.4 bindrcpp_0.2